rm(list = ls())
#working_folder<- "/Users/bgpopescu/Library/CloudStorage/Dropbox/covid_institutions/"
working_folder<-"//Mac/Dropbox-1/covid_paper_replication/"
setwd(working_folder)

#This is to obtain:
#-table_A7

library("readxl")
library("reshape")
library("stringi")
library("stringr")
library("dplyr")
library("tidyr")
library("stargazer")
library("corrplot")
library("Hmisc")
library("ggplot2")
library("glue")
library('gridExtra')
library('grid')
library('lfe')
library('broom')
library("plm")
library("multiwayvcov")
library("lmtest")
library("psych")
library("ggrepel")
library("lemon")
library("sandwich")
library("spdep")
library("xtable")
library('gsynth')
library("pracma")
library("conleyreg")
library("vars")
library("rlist")
library("fastDummies")
library("dplyr")
library("berryFunctions")
library("mediation")
library("sp")
library("pgirmess")
library("fwildclusterboot")
library("geojsonsf")
library("readr")


#Functions

dlshape=function(shploc, shpfile) {
  temp=tempfile()
  download.file(shploc, temp)
  unzip(temp)
  fp <- file.path(temp)
  obj_shp<-rgdal::readOGR(".",shpfile)
  return(obj_shp)}

mean_sd_fun<-function(lm_model) {
  list_empty <-c()
  number_mean<-mean(data.frame(lm_model$model)$mean_mobil, na.rm=T)
  number_sd<-sd(data.frame(lm_model$model)$mean_mobil, na.rm=T)
  list_empty[1]<-number_mean
  list_empty[2]<-number_sd
  names(list_empty)<-c("mean", "sd")
  return(list_empty)
}


week_obs_count<-function(lm_model) {
  list_empty <-c()
  week_no<-length(grep(x = colnames(data.frame(lm_model$model)), pattern = "^week"))
  #Adding 2 because there should be 53 weeks and because I am taking 2 weeks out as in STATA 2 weeks drop
  week_no<-week_no+2
  total_obs<-nobs(lm_model)
  #counties<-ncol(lm_model$model%>% dplyr:: select(starts_with("AGS_")))
  counties<-total_obs/week_no
  list_empty[1]<-week_no
  list_empty[2]<-total_obs
  list_empty[3]<-counties
  names(list_empty)<-c("week_no", "total_obs", "counties")
  return(list_empty)
}


#############################
#Function for calculating SE#
#############################

#Step1. Turning some variables to numeric
data_cov_mob = geojson_sf("./data/data_cov_mob.geojson")
data_cov_mob$pct_manufacturing<-as.numeric(data_cov_mob$pct_manufacturing)

data_cov_mob$kul<-as.numeric(data_cov_mob$kul)
data_cov_mob$nat<-as.numeric(data_cov_mob$nat)
data_cov_mob$frei<-as.numeric(data_cov_mob$frei)
data_cov_mob$soz<-as.numeric(data_cov_mob$soz)
data_cov_mob$pol<-as.numeric(data_cov_mob$pol)
data_cov_mob$inter<-as.numeric(data_cov_mob$inter)


data_cov_mob$gspo<-as.numeric(data_cov_mob$gspo)
data_cov_mob$gkul<-as.numeric(data_cov_mob$gkul)
data_cov_mob$gnat<-as.numeric(data_cov_mob$gnat)

data_cov_mob$gsoz<-as.numeric(data_cov_mob$gsoz)
data_cov_mob$gpol<-as.numeric(data_cov_mob$gpol)
data_cov_mob$ginter<-as.numeric(data_cov_mob$ginter)

#Bridging Networks
data_cov_mob$bridging_total<-data_cov_mob$gspo+ data_cov_mob$gkul+ data_cov_mob$gnat
data_cov_mob$bridging_nokul_total<-data_cov_mob$gspo+ data_cov_mob$gnat

data_cov_mob$bridging_tot_pop<-data_cov_mob$bridging_total/data_cov_mob$pop_total*1000
data_cov_mob$bridging_nokul_tot_pop<-data_cov_mob$bridging_nokul_total/data_cov_mob$pop_total*1000

mean_brid<-summary(data_cov_mob$bridging_tot_pop)[4]
mean_brid_nokul<-summary(data_cov_mob$bridging_nokul_tot_pop)[4]

data_cov_mob$bridging_tot_pop_dummy<-ifelse(data_cov_mob$bridging_tot_pop>mean_brid, 1, 0)
data_cov_mob$bridging_nokul_tot_pop_dummy<-ifelse(data_cov_mob$bridging_nokul_tot_pop>mean_brid_nokul, 1, 0)



#Bonding Networks
data_cov_mob$bonding_total<-data_cov_mob$gsoz+ data_cov_mob$gpol+ data_cov_mob$ginter
data_cov_mob$bonding_kul_total<-data_cov_mob$gsoz+ data_cov_mob$gpol+ data_cov_mob$ginter+ data_cov_mob$gkul

data_cov_mob$bonding_tot_pop<-data_cov_mob$bonding_total/data_cov_mob$pop_total*1000
data_cov_mob$bonding_kul_tot_pop<-data_cov_mob$bonding_kul_total/data_cov_mob$pop_total*1000

summary(data_cov_mob$bonding_tot_pop)
mean_bond<-summary(data_cov_mob$bonding_tot_pop)[4]
mean_bond_kul<-summary(data_cov_mob$bonding_kul_tot_pop)[4]

data_cov_mob$bonding_tot_pop_dummy<-ifelse(data_cov_mob$bonding_tot_pop>mean_bond, 1, 0)
data_cov_mob$bonding_kul_tot_pop_dummy<-ifelse(data_cov_mob$bonding_kul_tot_pop>mean_bond, 1, 0)


#All Networks
data_cov_mob$gvereine<-as.numeric(data_cov_mob$gvereine)
data_cov_mob$vereine_tot_pop<-data_cov_mob$gvereine/data_cov_mob$pop_total*1000
mean_ver<-summary(data_cov_mob$vereine_tot_pop)[4]
data_cov_mob$vereine_tot_pop_dummy<-ifelse(data_cov_mob$vereine_tot_pop>mean_ver, 1, 0)


data_cov_mob$gvereine<-as.numeric(data_cov_mob$gvereine)
data_cov_mob$vereine_tot_pop<-data_cov_mob$gvereine/data_cov_mob$pop_total*1000
mean_ver<-summary(data_cov_mob$vereine_tot_pop)[4]
data_cov_mob$vereine_tot_pop_dummy<-ifelse(data_cov_mob$vereine_tot_pop>mean_ver, 1, 0)


#Right Wing
data_cov_mob$pct_afd<-as.numeric(data_cov_mob$pct_afd)
mean_pct_afd<-summary(data_cov_mob$pct_afd)[3]
data_cov_mob$pct_afd_dummy<-ifelse(data_cov_mob$pct_afd>mean_pct_afd, 1, 0)

#Right Wing Change
data_cov_mob$pct_change_afd<-as.numeric(data_cov_mob$pct_change_afd)
mean_pct_change_afd<-summary(data_cov_mob$pct_change_afd)[4]
data_cov_mob$pct_change_afd_dummy<-ifelse(data_cov_mob$pct_afd>mean_pct_change_afd, 1, 0)


data_cov_mob_week<-subset(data_cov_mob, week=="01")
data_cov_mob_week<-subset(data_cov_mob_week, select = c("GEN", 
                                                        "vereine_tot_pop",
                                                        "bridging_tot_pop",
                                                        "bridging_nokul_tot_pop",
                                                        "bonding_kul_tot_pop",
                                                        "bonding_tot_pop",
                                                        "pct_afd_dummy",
                                                        "pct_change_afd_dummy",
                                                        "bridging_nokul_tot_pop"))
summary(data_cov_mob$gender_ratio)
data<-data_cov_mob
data$week_new<-as.numeric(data$week)+1
data$week_new<-paste0("0", data$week_new)  
data$week_new2<-ifelse(nchar(data$week_new)==3, str_remove(data$week_new, "^0+"), data$week_new)
data$week<-data$week_new2
data<-subset(data, select = -c(week_new, week_new2))
data$post_treat<-ifelse(as.numeric(data$week)>10, 1,0)

#data$bridging_nokul_tot_pop
data$post_treat_vereine_tot_pop_dummy<-data$post_treat*data$vereine_tot_pop_dummy

data$post_treat_bridging_tot_pop_dummy<-data$post_treat*data$bridging_tot_pop_dummy
data$post_treat_bridging_nokul_tot_pop_dummy<-data$post_treat*data$bridging_nokul_tot_pop_dummy

data$post_treat_bridging_tot_pop<-data$post_treat*data$bridging_tot_pop

data$post_treat_bonding_tot_pop_dummy<-data$post_treat*data$bonding_tot_pop_dummy
data$post_treat_bonding_kul_tot_pop_dummy<-data$post_treat*data$bonding_kul_tot_pop_dummy


data$post_treat_pct_afd_dummy<-data$post_treat*data$pct_afd_dummy
data$post_treat_pct_change_afd_dummy<-data$post_treat*data$pct_change_afd_dummy
data$post_treat_pct_change_afd<-data$post_treat*data$pct_change_afd


#Save to separate file
nuts<-subset(data_cov_mob, week=="00")
nuts2<-subset(nuts, select = c("AGS"))
nuts2$AGS<-as.numeric(nuts2$AGS<-nuts2$AGS)

#For Server too
# Create dummy variables:
dataf <- dummy_cols(data, select_columns = 'AGS')
dataw <- dummy_cols(dataf, select_columns = 'week')
datas <- dummy_cols(dataw, select_columns = 'SN_L')


new<-dataw %>% dplyr::select(starts_with("AGS_"))
new2<-dataw %>% dplyr::select(starts_with("week_"))
new3<-datas %>% dplyr::select(starts_with("SN_L"))


list_dummies<-names(new)
#I am removing the counties which get dropped in Stata

#Re-enable this later if necessary
list_dummies<-list_dummies[list_dummies!= "AGS_0"]
list_dummies<-list_dummies[list_dummies!= "AGS_5316"]
list_dummies<-list_dummies[list_dummies!= "AGS_8121"]
list_dummies<-list_dummies[list_dummies!= "AGS_8221"]
list_dummies<-list_dummies[list_dummies!= "AGS_8231"]
list_dummies<-list_dummies[list_dummies!= "AGS_8311"]
list_dummies<-list_dummies[list_dummies!= "AGS_8421"]
list_dummies<-list_dummies[list_dummies!= "AGS_9461"]
list_dummies<-list_dummies[list_dummies!= "AGS_13004"]
list_dummies<-list_dummies[list_dummies!= "AGS_16052"]
list_dummies<-list_dummies[list_dummies!= "AGS_16053"]
list_dummies<-list_dummies[list_dummies!= "AGS_16054"]
list_dummies<-list_dummies[list_dummies!= "AGS_16055"]
list_dummies<-list_dummies[list_dummies!= "AGS_16061"]
list_dummies<-list_dummies[list_dummies!= "AGS_16077"]

list_week_dummies<-names(new2)
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new2"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_numeric"]
list_week_dummies_now1<-list_week_dummies[list_week_dummies!= "week_01"]

#Re-enable this later if necessary
list_week_dummies2<-list_week_dummies[list_week_dummies!= "week_10"]
list_week_dummies2<-list_week_dummies2[list_week_dummies2!= "week_53"]

list_land_dummies<-names(new3)
#Re-enable this later if necessary
list_land_dummies_full<-list_land_dummies[list_land_dummies!= "SN_L"]
list_land_dummies<-list_land_dummies_full[list_land_dummies_full!= "SN_L_01"]
temp_lag<-10


#Creating the interaction terms for bonding
for (i in list_land_dummies){
  for (j in list_week_dummies_now1)
    datas[[paste("inter_", i, "_", j, sep="")]]<-datas[[i]]*datas[[j]]
}

new<-datas %>% dplyr::select(starts_with("inter_"))
list_land_dummies_star<-names(new)


#BRIDGING
basic_brid<-c("bridging_tot_pop_dummy", "post_treat", "post_treat_bridging_tot_pop_dummy", "lag_covid_pc")
extended_brid <- c("bridging_tot_pop_dummy", "post_treat", "post_treat_bridging_tot_pop_dummy", "lag_covid_pc",
                   "pct_turnout", "log_pop_total", "log_gdp_per_cap", "pop_den",
                   "pct_college", "pct_pop_above_60", "pct_pop_under_35", "pct_service_sec_narrow", "pct_manufacturing",
                   "east_ger", "gender_ratio", "pct_students")

basic_brid_nokul<-c("bridging_nokul_tot_pop_dummy", "post_treat", "post_treat_bridging_nokul_tot_pop_dummy", "lag_covid_pc")
extended_brid_nokul <- c("bridging_nokul_tot_pop_dummy", "post_treat", "post_treat_bridging_nokul_tot_pop_dummy", "lag_covid_pc",
                   "pct_turnout", "log_pop_total", "log_gdp_per_cap", "pop_den",
                   "pct_college", "pct_pop_above_60", "pct_pop_under_35", "pct_service_sec_narrow", "pct_manufacturing",
                   "east_ger", "gender_ratio", "pct_students")


#BONDING
basic_bond<-c("bonding_tot_pop_dummy", "post_treat", "post_treat_bonding_tot_pop_dummy", "lag_covid_pc")
extended_bond <- c("bonding_tot_pop_dummy", "post_treat", "post_treat_bonding_tot_pop_dummy", "lag_covid_pc",
                   "pct_turnout", "log_pop_total", "log_gdp_per_cap", "pop_den",
                   "pct_college", "pct_pop_above_60", "pct_pop_under_35", "pct_service_sec_narrow", "pct_manufacturing",
                   "east_ger", "gender_ratio", "pct_students")

basic_bond_kul<-c("bonding_kul_tot_pop_dummy", "post_treat", "post_treat_bonding_kul_tot_pop_dummy", "lag_covid_pc")
extended_bond_kul <- c("bonding_kul_tot_pop_dummy", "post_treat", "post_treat_bonding_kul_tot_pop_dummy", "lag_covid_pc",
                         "pct_turnout", "log_pop_total", "log_gdp_per_cap", "pop_den",
                         "pct_college", "pct_pop_above_60", "pct_pop_under_35", "pct_service_sec_narrow", "pct_manufacturing",
                         "east_ger", "gender_ratio", "pct_students")




basic_and_dummies_brid<-list.append(basic_brid, list_dummies, list_week_dummies2)
basic_and_dummies_int_brid<-list.append(basic_brid, list_dummies, list_week_dummies2, list_land_dummies_star)

basic_and_dummies_brid_nokul<-list.append(basic_brid_nokul, list_dummies, list_week_dummies2)
basic_and_dummies_int_brid_nokul<-list.append(basic_brid_nokul, list_dummies, list_week_dummies2, list_land_dummies_star)


extended_and_week_dummies_brid<-list.append(extended_brid, list_land_dummies, list_week_dummies2)
extended_and_week_dummies_int_brid<-list.append(extended_brid, list_land_dummies, list_week_dummies2, list_land_dummies_star)

extended_and_week_dummies_brid_nokul<-list.append(extended_brid_nokul, list_land_dummies, list_week_dummies2)
extended_and_week_dummies_int_brid_nokul<-list.append(extended_brid_nokul, list_land_dummies, list_week_dummies2, list_land_dummies_star)



basic_and_dummies_bond<-list.append(basic_bond, list_dummies, list_week_dummies2)
basic_and_dummies_int_bond<-list.append(basic_bond, list_dummies, list_week_dummies2, list_land_dummies_star)

basic_and_dummies_bond_kul<-list.append(basic_bond_kul, list_dummies, list_week_dummies2)
basic_and_dummies_int_bond_kul<-list.append(basic_bond_kul, list_dummies, list_week_dummies2, list_land_dummies_star)


extended_and_week_dummies_bond<-list.append(extended_bond, list_land_dummies, list_week_dummies2)
extended_and_week_dummies_int_bond<-list.append(extended_bond, list_land_dummies, list_week_dummies2, list_land_dummies_star)

extended_and_week_dummies_bond_kul<-list.append(extended_bond_kul, list_land_dummies, list_week_dummies2)
extended_and_week_dummies_int_bond_kul<-list.append(extended_bond_kul, list_land_dummies, list_week_dummies2, list_land_dummies_star)


spec_basic_brid<-as.formula(paste("mean_mobil~",
                                  paste(basic_and_dummies_brid, collapse= "+")))
spec_basic_int_brid<-as.formula(paste("mean_mobil~",
                                      paste(basic_and_dummies_int_brid, collapse= "+")))

spec_extended_brid<-as.formula(paste("mean_mobil~",
                                     paste(extended_and_week_dummies_brid, collapse= "+")))
spec_extended_int_brid<-as.formula(paste("mean_mobil~",
                                         paste(extended_and_week_dummies_int_brid, collapse= "+")))

spec_basic_brid_nokul<-as.formula(paste("mean_mobil~",
                                  paste(basic_and_dummies_brid_nokul, collapse= "+")))
spec_basic_int_brid_nokul<-as.formula(paste("mean_mobil~",
                                      paste(basic_and_dummies_int_brid_nokul, collapse= "+")))

spec_extended_brid_nokul<-as.formula(paste("mean_mobil~",
                                     paste(extended_and_week_dummies_brid_nokul, collapse= "+")))
spec_extended_int_brid_nokul<-as.formula(paste("mean_mobil~",
                                         paste(extended_and_week_dummies_int_brid_nokul, collapse= "+")))



spec_basic_bond<-as.formula(paste("mean_mobil~",
                                  paste(basic_and_dummies_bond, collapse= "+")))
spec_basic_int_bond<-as.formula(paste("mean_mobil~",
                                      paste(basic_and_dummies_int_bond, collapse= "+")))

spec_extended_bond<-as.formula(paste("mean_mobil~",
                                     paste(extended_and_week_dummies_bond, collapse= "+")))
spec_extended_int_bond<-as.formula(paste("mean_mobil~",
                                         paste(extended_and_week_dummies_int_bond, collapse= "+")))

spec_basic_bond_kul<-as.formula(paste("mean_mobil~",
                                      paste(basic_and_dummies_bond_kul, collapse= "+")))
spec_basic_int_bond_kul<-as.formula(paste("mean_mobil~",
                                          paste(basic_and_dummies_int_bond_kul, collapse= "+")))

spec_extended_bond_kul<-as.formula(paste("mean_mobil~",
                                         paste(extended_and_week_dummies_bond_kul, collapse= "+")))
spec_extended_int_bond_kul<-as.formula(paste("mean_mobil~",
                                             paste(extended_and_week_dummies_int_bond_kul, collapse= "+")))


################
#Basic Bridging#
################
mo1_brid<-lm(spec_basic_int_brid, data=datas)
mo1_brid_df<-broom::tidy(mo1_brid)
mo1_brid_df$rse<-coeftest(mo1_brid, vcov = cluster.vcov(mo1_brid, datas$AGS, force_posdef=T))[, 2]
mo1_brid_df$pv_rse<-coeftest(mo1_brid, vcov = cluster.vcov(mo1_brid, datas$AGS, force_posdef=T))[, 4]
mo1_brid_df2<-mo1_brid_df %>% filter(term %in% basic_brid)

basic_and_dummies_int_brid_x<-list.append(basic_brid, list_land_dummies_star)
spec_basic_int_brid_x<-as.formula(paste("mean_mobil~",
                                        paste(basic_and_dummies_int_brid_x, collapse= "+"), "| AGS + week"))


x3<-labels(terms(spec_basic_int_brid))
datas$time<-as.numeric(datas$week)

datas3<-subset(datas, select = list.append(x3, "bridging_total", "mean_mobil", "POINT_X", "POINT_Y", "time", "AGS"))
datas3<-datas3[complete.cases(datas3), ]
datas3_oneweek<-subset(datas3, week_11==1)
datas3_oneweek$AGS<-as.numeric(datas3_oneweek$AGS)

#Conley SE

#Deciding on the temporal lag - value with the lowest AIC
x<-VARselect(datas3$mean_mobil)
min(x$criteria["AIC(n)",])
temp_lag<-names(x$criteria["AIC(n)",])[x$criteria["AIC(n)",]==min(x$criteria["AIC(n)",])]
temp_lag<-as.numeric(temp_lag)


#Merging with spatial dataframe
datas3_oneweek<-st_drop_geometry(datas3_oneweek)
nuts_merged<-left_join(nuts2, datas3_oneweek, by = c("AGS"="AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$bridging_total))

nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

coo<-subset(nuts_merged3, select = c("POINT_X_m", "POINT_Y_m"))
coo<-st_drop_geometry(coo)


#Calculating the correlogram
pgi_cor <- correlog(coords=coo, z=nuts_merged3$bridging_total, method="Moran", 
                    nbclass=10)
pgi_cor_df<-data.frame(pgi_cor)



#Calculating the distance the minimises Moran's I stat
dist_cut_brid<-round(pgi_cor_df$dist.class[abs(pgi_cor_df$coef)== min(abs(pgi_cor_df$coef))]/100000,0)

dm_bri <- dist_mat(datas3, dist_cutoff = dist_cut_brid, 
                   lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time",
                   dist_round = TRUE)
res<-conleyreg(spec_basic_int_brid, data=datas3, dist_cutoff=dist_cut_brid, unit = "AGS", time = "time", 
               dist_mat=dm_bri, lag_cutoff=temp_lag)

res_basic_int_brid<-broom::tidy(res)
res_basic_int_brid_df2<-res_basic_int_brid %>% filter(term %in% basic_brid)
res_basic_int_brid_df2<-subset(res_basic_int_brid_df2, select = c("term", "std.error", 'p.value'))
names(res_basic_int_brid_df2)<-c("term", "con_std_error", 'con_p_value')

model3<-left_join(mo1_brid_df2, res_basic_int_brid_df2, by = c("term" = "term"))
dir.create(file.path(working_folder, "./paper/tables/tableA6/"), showWarnings = FALSE)
write.csv(model3,"./paper/tables/tableA6/model3_brid.csv", row.names = FALSE)


# *Model4
# *NO County FE
mo8_brid<-lm(spec_extended_int_brid, data=datas)
summary(mo8_brid)

mo8_brid_df<-broom::tidy(mo8_brid)
mo8_brid_df$rse<-coeftest(mo8_brid, vcov = cluster.vcov(mo8_brid, datas$AGS, force_posdef=T))[, 2]
mo8_brid_df$pv_rse<-coeftest(mo8_brid, vcov = cluster.vcov(mo8_brid, datas$AGS, force_posdef=T))[, 4]
mo8_brid_df2<-mo8_brid_df %>% filter(term %in% extended_brid)


x3<-labels(terms(spec_extended_int_brid))
datas$time<-as.numeric(datas$week)

datas3<-subset(datas, select = list.append(x3, "mean_mobil", "bridging_total", "POINT_X", "POINT_Y", "time", "AGS"))
datas3<-datas3[complete.cases(datas3), ]

dm_bri <- dist_mat(datas3, dist_cutoff = dist_cut_brid, 
                   lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time",
                   dist_round = TRUE)
res<-conleyreg(spec_extended_int_brid, data=datas3, dist_cutoff=dist_cut_brid, unit = "AGS", time = "time", 
               dist_mat=dm_bri, lag_cutoff=temp_lag)

res_extended_int_brid_df<-broom::tidy(res)
res_extended_int_brid_df2<-res_extended_int_brid_df %>% filter(term %in% extended_brid)
res_extended_int_brid_df2<-subset(res_extended_int_brid_df, select = c("term", "std.error", 'p.value'))
names(res_extended_int_brid_df2)<-c("term", "con_std_error", 'con_p_value')

model4<-left_join(mo8_brid_df2, res_extended_int_brid_df2, by = c("term" = "term"))
write.csv(model4,"./paper/tables/tableA6/model4_brid.csv", row.names = FALSE)

################
#Basic Bonding#
################
mo1_bond<-lm(spec_basic_int_bond, data=datas)
mo1_bond_df<-broom::tidy(mo1_bond)
mo1_bond_df$rse<-coeftest(mo1_bond, vcov = cluster.vcov(mo1_bond, datas$AGS, force_posdef=T))[, 2]
mo1_bond_df$pv_rse<-coeftest(mo1_bond, vcov = cluster.vcov(mo1_bond, datas$AGS, force_posdef=T))[, 4]
mo1_bond_df2<-mo1_bond_df %>% filter(term %in% basic_bond)

basic_and_dummies_int_bond_x<-list.append(basic_bond, list_land_dummies_star)
spec_basic_int_bond_x<-as.formula(paste("mean_mobil~",
                                        paste(basic_and_dummies_int_bond_x, collapse= "+"), "| AGS + week"))

x3<-labels(terms(spec_basic_int_bond))
datas$time<-as.numeric(datas$week)

datas3<-subset(datas, select = list.append(x3, "bonding_total", "mean_mobil", "POINT_X", "POINT_Y", "time", "AGS"))
datas3<-datas3[complete.cases(datas3), ]
datas3_oneweek<-subset(datas3, week_11==1)
datas3_oneweek$AGS<-as.numeric(datas3_oneweek$AGS)

#Conley SE

#Deciding on the temporal lag - value with the lowest AIC
x<-VARselect(datas3$mean_mobil)
min(x$criteria["AIC(n)",])
temp_lag<-names(x$criteria["AIC(n)",])[x$criteria["AIC(n)",]==min(x$criteria["AIC(n)",])]
temp_lag<-as.numeric(temp_lag)


#Merging with spatial dataframe
datas3_oneweek<-st_drop_geometry(datas3_oneweek)
nuts_merged<-left_join(nuts2, datas3_oneweek, by = c("AGS"="AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$bonding_total))

nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

coo<-subset(nuts_merged3, select = c("POINT_X_m", "POINT_Y_m"))
coo<-st_drop_geometry(coo)


#Calculating the correlogram
pgi_cor <- correlog(coords=coo, z=nuts_merged3$bonding_total, method="Moran", 
                    nbclass=10)
pgi_cor_df<-data.frame(pgi_cor)



#Calculating the distance the minimises Moran's I stat
dist_cut_bond<-round(pgi_cor_df$dist.class[abs(pgi_cor_df$coef)== min(abs(pgi_cor_df$coef))]/100000,0)

dm_bri <- dist_mat(datas3, dist_cutoff = dist_cut_bond, 
                   lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time",
                   dist_round = TRUE)
res<-conleyreg(spec_basic_int_bond, data=datas3, dist_cutoff=dist_cut_bond, unit = "AGS", time = "time", 
               dist_mat=dm_bri, lag_cutoff=temp_lag)

res_basic_int_bond<-broom::tidy(res)
res_basic_int_bond_df2<-res_basic_int_bond %>% filter(term %in% basic_bond)
res_basic_int_bond_df2<-subset(res_basic_int_bond_df2, select = c("term", "std.error", 'p.value'))
names(res_basic_int_bond_df2)<-c("term", "con_std_error", 'con_p_value')

model3<-left_join(mo1_bond_df2, res_basic_int_bond_df2, by = c("term" = "term"))
write.csv(model3,"./paper/tables/tableA6/model3_bond.csv", row.names = FALSE)


# *Model4
# *NO County FE
mo8_bond<-lm(spec_extended_int_bond, data=datas)
summary(mo8_bond)

mo8_bond_df<-broom::tidy(mo8_bond)
mo8_bond_df$rse<-coeftest(mo8_bond, vcov = cluster.vcov(mo8_bond, datas$AGS, force_posdef=T))[, 2]
mo8_bond_df$pv_rse<-coeftest(mo8_bond, vcov = cluster.vcov(mo8_bond, datas$AGS, force_posdef=T))[, 4]
mo8_bond_df2<-mo8_bond_df %>% filter(term %in% extended_bond)


x3<-labels(terms(spec_extended_int_bond))
datas$time<-as.numeric(datas$week)

datas3<-subset(datas, select = list.append(x3, "mean_mobil", "bonding_total", "POINT_X", "POINT_Y", "time", "AGS"))
datas3<-datas3[complete.cases(datas3), ]


dm_bri <- dist_mat(datas3, dist_cutoff = dist_cut_bond, 
                   lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time",
                   dist_round = TRUE)
res<-conleyreg(spec_extended_int_bond, data=datas3, dist_cutoff=dist_cut_bond, unit = "AGS", time = "time", 
               dist_mat=dm_bri, lag_cutoff=temp_lag)

res_extended_int_bond_df<-broom::tidy(res)
res_extended_int_bond_df2<-res_extended_int_bond_df %>% filter(term %in% extended_bond)
res_extended_int_bond_df2<-subset(res_extended_int_bond_df, select = c("term", "std.error", 'p.value'))
names(res_extended_int_bond_df2)<-c("term", "con_std_error", 'con_p_value')

model4<-left_join(mo8_bond_df2, res_extended_int_bond_df2, by = c("term" = "term"))
write.csv(model4,"./paper/tables/tableA6/model4_bond.csv", row.names = FALSE)



#######################
#Basic Bridging NO KUL#
#######################
mo1_brid_nokul<-lm(spec_basic_int_brid_nokul, data=datas)
mo1_brid_nokul_df<-broom::tidy(mo1_brid_nokul)
mo1_brid_nokul_df$rse<-coeftest(mo1_brid_nokul, vcov = cluster.vcov(mo1_brid_nokul, datas$AGS, force_posdef=T))[, 2]
mo1_brid_nokul_df$pv_rse<-coeftest(mo1_brid_nokul, vcov = cluster.vcov(mo1_brid_nokul, datas$AGS, force_posdef=T))[, 4]
mo1_brid_nokul_df2<-mo1_brid_nokul_df %>% filter(term %in% basic_brid_nokul)


x3<-labels(terms(spec_basic_int_brid_nokul))
datas$time<-as.numeric(datas$week)

datas3<-subset(datas, select = list.append(x3, "bridging_nokul_total", "mean_mobil", "POINT_X", "POINT_Y", "time", "AGS"))
datas3<-datas3[complete.cases(datas3), ]
datas3_oneweek<-subset(datas3, week_11==1)
datas3_oneweek$AGS<-as.numeric(datas3_oneweek$AGS)

#Conley SE

#Deciding on the temporal lag - value with the lowest AIC
x<-VARselect(datas3$mean_mobil)
min(x$criteria["AIC(n)",])
temp_lag<-names(x$criteria["AIC(n)",])[x$criteria["AIC(n)",]==min(x$criteria["AIC(n)",])]
temp_lag<-as.numeric(temp_lag)


#Merging with spatial dataframe
datas3_oneweek<-st_drop_geometry(datas3_oneweek)
datas3_oneweek$bridging_nokul_total
nuts_merged<-left_join(nuts2, datas3_oneweek, by = c("AGS"="AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$bridging_nokul_total))

nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

coo<-subset(nuts_merged3, select = c("POINT_X_m", "POINT_Y_m"))
coo<-st_drop_geometry(coo)


#Calculating the correlogram
pgi_cor <- correlog(coords=coo, z=nuts_merged3$bridging_nokul_total, method="Moran", 
                    nbclass=10)
pgi_cor_df<-data.frame(pgi_cor)



#Calculating the distance the minimises Moran's I stat
dist_cut_brid<-round(pgi_cor_df$dist.class[abs(pgi_cor_df$coef)== min(abs(pgi_cor_df$coef))]/100000,0)

dm_bri <- dist_mat(datas3, dist_cutoff = dist_cut_brid, 
                   lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time",
                   dist_round = TRUE)
res<-conleyreg(spec_basic_int_brid_nokul, data=datas3, dist_cutoff=dist_cut_brid, unit = "AGS", time = "time", 
               dist_mat=dm_bri, lag_cutoff=temp_lag)

res_basic_int_brid_nokul<-broom::tidy(res)
res_basic_int_brid_nokul_df2<-res_basic_int_brid_nokul %>% filter(term %in% basic_brid_nokul)
res_basic_int_brid_nokul_df2<-subset(res_basic_int_brid_nokul_df2, select = c("term", "std.error", 'p.value'))
names(res_basic_int_brid_nokul_df2)<-c("term", "con_std_error", 'con_p_value')


model3_nokul<-left_join(mo1_brid_nokul_df2, res_basic_int_brid_nokul_df2, by = c("term" = "term"))
write.csv(model3_nokul,"./paper/tables/tableA6/model3_brid_nokul.csv", row.names = FALSE)


# *Model4
# *NO County FE
mo8_brid_nokul<-lm(spec_extended_int_brid_nokul, data=datas)
summary(mo8_brid_nokul)

mo8_brid_nokul_df<-broom::tidy(mo8_brid_nokul)
mo8_brid_nokul_df$rse<-coeftest(mo8_brid_nokul, vcov = cluster.vcov(mo8_brid_nokul, datas$AGS, force_posdef=T))[, 2]
mo8_brid_nokul_df$pv_rse<-coeftest(mo8_brid_nokul, vcov = cluster.vcov(mo8_brid_nokul, datas$AGS, force_posdef=T))[, 4]
mo8_brid_nokul_df2<-mo8_brid_nokul_df %>% filter(term %in% extended_brid_nokul)


x3<-labels(terms(spec_extended_int_brid_nokul))
datas$time<-as.numeric(datas$week)

datas3<-subset(datas, select = list.append(x3, "mean_mobil", "bridging_nokul_total", "POINT_X", "POINT_Y", "time", "AGS"))
datas3<-datas3[complete.cases(datas3), ]

dm_bri <- dist_mat(datas3, dist_cutoff = dist_cut_brid, 
                   lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time",
                   dist_round = TRUE)

datas3$bridging_nokul_tot_pop_dummy
datas3$post_treat
datas3$post_treat_bridging_nokul_tot_pop_dummy
datas3$lag_covid_pc

res<-conleyreg(spec_extended_int_brid_nokul, data=datas3, dist_cutoff=dist_cut_brid, unit = "AGS", time = "time", 
               dist_mat=dm_bri, lag_cutoff=temp_lag)

res_extended_int_brid_df<-broom::tidy(res)
res_extended_int_brid_df2<-res_extended_int_brid_df %>% filter(term %in% extended_brid_nokul)
res_extended_int_brid_df2<-subset(res_extended_int_brid_df, select = c("term", "std.error", 'p.value'))
names(res_extended_int_brid_df2)<-c("term", "con_std_error", 'con_p_value')

model4<-left_join(mo8_brid_nokul_df2, res_extended_int_brid_df2, by = c("term" = "term"))
write.csv(model4,"./paper/tables/tableA6/model4_brid_nokul.csv", row.names = FALSE)



#######################
#Basic Bonding KUL#
#######################
mo1_bond_kul<-lm(spec_basic_int_bond_kul, data=datas)
mo1_bond_kul_df<-broom::tidy(mo1_bond_kul)
mo1_bond_kul_df$rse<-coeftest(mo1_bond_kul, vcov = cluster.vcov(mo1_bond_kul, datas$AGS, force_posdef=T))[, 2]
mo1_bond_kul_df$pv_rse<-coeftest(mo1_bond_kul, vcov = cluster.vcov(mo1_bond_kul, datas$AGS, force_posdef=T))[, 4]
mo1_bond_kul_df2<-mo1_bond_kul_df %>% filter(term %in% basic_bond_kul)


x3<-labels(terms(spec_basic_int_bond_kul))
datas$time<-as.numeric(datas$week)

datas3<-subset(datas, select = list.append(x3, "bridging_nokul_total", "mean_mobil", "POINT_X", "POINT_Y", "time", "AGS"))
datas3<-datas3[complete.cases(datas3), ]
datas3_oneweek<-subset(datas3, week_11==1)
datas3_oneweek$AGS<-as.numeric(datas3_oneweek$AGS)

#Conley SE

#Deciding on the temporal lag - value with the lowest AIC
x<-VARselect(datas3$mean_mobil)
min(x$criteria["AIC(n)",])
temp_lag<-names(x$criteria["AIC(n)",])[x$criteria["AIC(n)",]==min(x$criteria["AIC(n)",])]
temp_lag<-as.numeric(temp_lag)


#Merging with spatial dataframe
datas3_oneweek<-st_drop_geometry(datas3_oneweek)
datas3_oneweek$bridging_nokul_total
nuts_merged<-left_join(nuts2, datas3_oneweek, by = c("AGS"="AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$bridging_nokul_total))

nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

coo<-subset(nuts_merged3, select = c("POINT_X_m", "POINT_Y_m"))
coo<-st_drop_geometry(coo)


#Calculating the correlogram
pgi_cor <- correlog(coords=coo, z=nuts_merged3$bridging_nokul_total, method="Moran", 
                    nbclass=10)
pgi_cor_df<-data.frame(pgi_cor)



#Calculating the distance the minimises Moran's I stat
dist_cut_brid<-round(pgi_cor_df$dist.class[abs(pgi_cor_df$coef)== min(abs(pgi_cor_df$coef))]/100000,0)

dm_bri <- dist_mat(datas3, dist_cutoff = dist_cut_brid, 
                   lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time",
                   dist_round = TRUE)
res<-conleyreg(spec_basic_int_bond_kul, data=datas3, dist_cutoff=dist_cut_brid, unit = "AGS", time = "time", 
               dist_mat=dm_bri, lag_cutoff=temp_lag)

res_basic_int_bond_kul<-broom::tidy(res)
res_basic_int_bond_kul_df2<-res_basic_int_bond_kul %>% filter(term %in% basic_bond_kul)
res_basic_int_bond_kul_df2<-subset(res_basic_int_bond_kul_df2, select = c("term", "std.error", 'p.value'))
names(res_basic_int_bond_kul_df2)<-c("term", "con_std_error", 'con_p_value')


model3_nokul<-left_join(mo1_bond_kul_df2, res_basic_int_bond_kul_df2, by = c("term" = "term"))
write.csv(model3_nokul,"./paper/tables/tableA6/model3_bond_kul.csv", row.names = FALSE)


# *Model4
# *NO County FE
mo8_bond_kul<-lm(spec_extended_int_bond_kul, data=datas)
summary(mo8_bond_kul)

mo8_bond_kul_df<-broom::tidy(mo8_bond_kul)
mo8_bond_kul_df$rse<-coeftest(mo8_bond_kul, vcov = cluster.vcov(mo8_bond_kul, datas$AGS, force_posdef=T))[, 2]
mo8_bond_kul_df$pv_rse<-coeftest(mo8_bond_kul, vcov = cluster.vcov(mo8_bond_kul, datas$AGS, force_posdef=T))[, 4]
mo8_bond_kul_df2<-mo8_bond_kul_df %>% filter(term %in% extended_bond_kul)


x3<-labels(terms(spec_extended_int_bond_kul))
datas$time<-as.numeric(datas$week)

datas3<-subset(datas, select = list.append(x3, "mean_mobil", "bridging_nokul_total", "POINT_X", "POINT_Y", "time", "AGS"))
datas3<-datas3[complete.cases(datas3), ]


dm_bri <- dist_mat(datas3, dist_cutoff = dist_cut_brid, 
                   lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time",
                   dist_round = TRUE)

datas3$bridging_nokul_tot_pop_dummy
datas3$post_treat
datas3$post_treat_bridging_nokul_tot_pop_dummy
datas3$lag_covid_pc

res<-conleyreg(spec_extended_int_bond_kul, data=datas3, dist_cutoff=dist_cut_brid, unit = "AGS", time = "time", 
               dist_mat=dm_bri, lag_cutoff=temp_lag)

res_extended_int_brid_df<-broom::tidy(res)
res_extended_int_brid_df2<-res_extended_int_brid_df %>% filter(term %in% extended_bond_kul)
res_extended_int_brid_df2<-subset(res_extended_int_brid_df, select = c("term", "std.error", 'p.value'))
names(res_extended_int_brid_df2)<-c("term", "con_std_error", 'con_p_value')

model4<-left_join(mo8_bond_kul_df2, res_extended_int_brid_df2, by = c("term" = "term"))
write.csv(model4,"./paper/tables/tableA6/model4_bond_kul.csv", row.names = FALSE)


model1_csv <- as.data.frame(read_csv("paper/tables/tableA6/model3_brid.csv"))
model2_csv <- as.data.frame(read_csv("paper/tables/tableA6/model4_brid.csv"))
model3_csv <- as.data.frame(read_csv("paper/tables/tableA6/model3_brid_nokul.csv"))
model4_csv <- as.data.frame(read_csv("paper/tables/tableA6/model4_brid_nokul.csv"))
model5_csv <- as.data.frame(read_csv("paper/tables/tableA6/model3_bond.csv"))
model6_csv <- as.data.frame(read_csv("paper/tables/tableA6/model4_bond.csv"))
model7_csv <- as.data.frame(read_csv("paper/tables/tableA6/model3_bond_kul.csv"))
model8_csv <- as.data.frame(read_csv("paper/tables/tableA6/model4_bond_kul.csv"))


model1<-model1_csv
model2<-model2_csv
model3<-model3_csv
model4<-model4_csv
model5<-model5_csv
model6<-model6_csv
model7<-model7_csv
model8<-model8_csv

model1$order<-NA
model2$order<-NA
model3$order<-NA
model4$order<-NA
model5$order<-NA
model6$order<-NA
model7$order<-NA
model8$order<-NA



model1$order[model1$term=="bridging_tot_pop_dummy"]<-1
model2$order[model2$term=="bridging_tot_pop_dummy"]<-1
model3$order[model3$term=="bridging_nokul_tot_pop_dummy"]<-2
model4$order[model4$term=="bridging_nokul_tot_pop_dummy"]<-2
model5$order[model5$term=="bonding_tot_pop_dummy"]<-3
model6$order[model6$term=="bonding_tot_pop_dummy"]<-3
model7$order[model7$term=="bonding_kul_tot_pop_dummy"]<-4
model8$order[model8$term=="bonding_kul_tot_pop_dummy"]<-4

model1$order[model1$term=="post_treat"]<-5
model2$order[model2$term=="post_treat"]<-5
model3$order[model3$term=="post_treat"]<-5
model4$order[model4$term=="post_treat"]<-5
model5$order[model5$term=="post_treat"]<-5
model6$order[model6$term=="post_treat"]<-5
model7$order[model7$term=="post_treat"]<-5
model8$order[model8$term=="post_treat"]<-5



model1$order[model1$term=="post_treat_bridging_tot_pop_dummy"]<-6
model2$order[model2$term=="post_treat_bridging_tot_pop_dummy"]<-6
model3$order[model3$term=="post_treat_bridging_nokul_tot_pop_dummy"]<-7
model4$order[model4$term=="post_treat_bridging_nokul_tot_pop_dummy"]<-7
model5$order[model5$term=="post_treat_bonding_tot_pop_dummy"]<-8
model6$order[model6$term=="post_treat_bonding_tot_pop_dummy"]<-8
model7$order[model7$term=="post_treat_bonding_kul_tot_pop_dummy"]<-9
model8$order[model8$term=="post_treat_bonding_kul_tot_pop_dummy"]<-9



model2$order[model2$term=="pct_turnout"]<-10
model4$order[model4$term=="pct_turnout"]<-10
model6$order[model6$term=="pct_turnout"]<-10
model8$order[model8$term=="pct_turnout"]<-10


model2$order[model2$term=="log_pop_total"]<-11
model4$order[model4$term=="log_pop_total"]<-11
model6$order[model6$term=="log_pop_total"]<-11
model8$order[model8$term=="log_pop_total"]<-11


model2$order[model2$term=="log_gdp_per_cap"]<-12
model4$order[model4$term=="log_gdp_per_cap"]<-12
model6$order[model6$term=="log_gdp_per_cap"]<-12
model8$order[model8$term=="log_gdp_per_cap"]<-12


model2$order[model2$term=="pop_den"]<-13
model4$order[model4$term=="pop_den"]<-13
model6$order[model6$term=="pop_den"]<-13
model8$order[model8$term=="pop_den"]<-13


model2$order[model2$term=="pct_college"]<-14
model4$order[model4$term=="pct_college"]<-14
model6$order[model6$term=="pct_college"]<-14
model8$order[model8$term=="pct_college"]<-14


model2$order[model2$term=="pct_pop_above_60"]<-15
model4$order[model4$term=="pct_pop_above_60"]<-15
model6$order[model6$term=="pct_pop_above_60"]<-15
model8$order[model8$term=="pct_pop_above_60"]<-15


model2$order[model2$term=="pct_pop_under_35"]<-16
model4$order[model4$term=="pct_pop_under_35"]<-16
model6$order[model6$term=="pct_pop_under_35"]<-16
model8$order[model8$term=="pct_pop_under_35"]<-16

model2$order[model2$term=="pct_service_sec_narrow"]<-17
model4$order[model4$term=="pct_service_sec_narrow"]<-17
model6$order[model6$term=="pct_service_sec_narrow"]<-17
model8$order[model8$term=="pct_service_sec_narrow"]<-17


model2$order[model2$term=="pct_manufacturing"]<-18
model4$order[model4$term=="pct_manufacturing"]<-18
model6$order[model6$term=="pct_manufacturing"]<-18
model8$order[model8$term=="pct_manufacturing"]<-18

model2$order[model2$term=="east_ger"]<-19
model4$order[model4$term=="east_ger"]<-19
model6$order[model6$term=="east_ger"]<-19
model8$order[model8$term=="east_ger"]<-19


model1$order[model1$term=="lag_covid_pc"]<-20
model2$order[model2$term=="lag_covid_pc"]<-20
model6$order[model6$term=="lag_covid_pc"]<-20
model8$order[model8$term=="lag_covid_pc"]<-20

model2$order[model2$term=="gender_ratio"]<-21
model4$order[model4$term=="gender_ratio"]<-21
model6$order[model6$term=="gender_ratio"]<-21
model8$order[model8$term=="gender_ratio"]<-21

model2$order[model2$term=="pct_students"]<-22
model4$order[model4$term=="pct_students"]<-22
model6$order[model6$term=="pct_students"]<-22
model8$order[model8$term=="pct_students"]<-22


model1$estimate<-round(as.numeric(model1$estimate), 3)
model2$estimate<-round(as.numeric(model2$estimate), 3)
model3$estimate<-round(as.numeric(model3$estimate), 3)
model4$estimate<-round(as.numeric(model4$estimate), 3)
model5$estimate<-round(as.numeric(model5$estimate), 3)
model6$estimate<-round(as.numeric(model6$estimate), 3)
model7$estimate<-round(as.numeric(model7$estimate), 3)
model8$estimate<-round(as.numeric(model8$estimate), 3)

model1$std_er_star<-paste("(", round(as.numeric(model1$rse),3), ")",
                          ifelse(as.numeric(model1$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model1$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model1$pv_rse)<0.05,"*",
                                               ""))), sep = "")

model2$std_er_star<-paste("(", round(as.numeric(model2$rse),3), ")",
                          ifelse(as.numeric(model2$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model2$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model2$pv_rse)<0.05,"*",
                                               ""))), sep = "")
model3$std_er_star<-paste("(", round(as.numeric(model3$rse),3), ")",
                          ifelse(as.numeric(model3$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model3$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model3$pv_rse)<0.05,"*",
                                               ""))), sep = "")

model4$std_er_star<-paste("(", round(as.numeric(model4$rse),3), ")",
                          ifelse(as.numeric(model4$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model4$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model4$pv_rse)<0.05,"*",
                                               ""))), sep = "")

model5$std_er_star<-paste("(", round(as.numeric(model5$rse),3), ")",
                          ifelse(as.numeric(model5$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model5$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model5$pv_rse)<0.05,"*",
                                               ""))), sep = "")

model6$std_er_star<-paste("(", round(as.numeric(model6$rse),3), ")",
                          ifelse(as.numeric(model6$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model6$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model6$pv_rse)<0.05,"*",
                                               ""))), sep = "")

model7$std_er_star<-paste("(", round(as.numeric(model7$rse),3), ")",
                          ifelse(as.numeric(model7$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model7$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model7$pv_rse)<0.05,"*",
                                               ""))), sep = "")


model8$std_er_star<-paste("(", round(as.numeric(model8$rse),3), ")",
                          ifelse(as.numeric(model8$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model8$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model8$pv_rse)<0.05,"*",
                                               ""))), sep = "")


model1$con_std_er_star<-paste("[", round(as.numeric(model1$con_std_error),3), "]",
                              ifelse(as.numeric(model1$con_p_value)<0.001,"***",
                                     ifelse(as.numeric(model1$con_p_value<0.01),"**",
                                            ifelse(as.numeric(model1$con_p_value)<0.05,"*",
                                                   ""))), sep = "")


model2$con_std_er_star<-paste("[", round(as.numeric(model2$con_std_error),3), "]",
                              ifelse(as.numeric(model2$con_p_value)<0.001,"***",
                                     ifelse(as.numeric(model2$con_p_value<0.01),"**",
                                            ifelse(as.numeric(model2$con_p_value)<0.05,"*",
                                                   ""))), sep = "")

model3$con_std_er_star<-paste("[", round(as.numeric(model3$con_std_error),3), "]",
                              ifelse(as.numeric(model3$con_p_value)<0.001,"***",
                                     ifelse(as.numeric(model3$con_p_value<0.01),"**",
                                            ifelse(as.numeric(model3$con_p_value)<0.05,"*",
                                                   ""))), sep = "")

model4$con_std_er_star<-paste("[", round(as.numeric(model4$con_std_error),3), "]",
                              ifelse(as.numeric(model4$con_p_value)<0.001,"***",
                                     ifelse(as.numeric(model4$con_p_value<0.01),"**",
                                            ifelse(as.numeric(model4$con_p_value)<0.05,"*",
                                                   ""))), sep = "")


model5$con_std_er_star<-paste("[", round(as.numeric(model5$con_std_error),3), "]",
                              ifelse(as.numeric(model5$con_p_value)<0.001,"***",
                                     ifelse(as.numeric(model5$con_p_value<0.01),"**",
                                            ifelse(as.numeric(model5$con_p_value)<0.05,"*",
                                                   ""))), sep = "")


model6$con_std_er_star<-paste("[", round(as.numeric(model6$con_std_error),3), "]",
                              ifelse(as.numeric(model6$con_p_value)<0.001,"***",
                                     ifelse(as.numeric(model6$con_p_value<0.01),"**",
                                            ifelse(as.numeric(model6$con_p_value)<0.05,"*",
                                                   ""))), sep = "")

model7$con_std_er_star<-paste("[", round(as.numeric(model7$con_std_error),3), "]",
                              ifelse(as.numeric(model7$con_p_value)<0.001,"***",
                                     ifelse(as.numeric(model7$con_p_value<0.01),"**",
                                            ifelse(as.numeric(model7$con_p_value)<0.05,"*",
                                                   ""))), sep = "")

model8$con_std_er_star<-paste("[", round(as.numeric(model8$con_std_error),3), "]",
                              ifelse(as.numeric(model8$con_p_value)<0.001,"***",
                                     ifelse(as.numeric(model8$con_p_value<0.01),"**",
                                            ifelse(as.numeric(model8$con_p_value)<0.05,"*",
                                                   ""))), sep = "")


mode1_b<-subset(model1, select = c("order", "term", "estimate", "std_er_star", "con_std_er_star"))
means<-c(23, "mean_mob", round(mean_sd_fun(mo1_brid)[1], 3), "", "", "")
sds<-c(24, "sd_mob", round(mean_sd_fun(mo1_brid)[2], 3), "", "", "")
no_cross_sections<-c(25, "cross_sections", round(week_obs_count(mo1_brid)[3], 3), "", "", "")
no_time_units<-c(26, "time_units", round(week_obs_count(mo1_brid)[1], 3), "", "", "")
county_fe<-c(27, "county_fe", "Yes", "", "", "")
week_fe<-c(28, "week_fe", "Yes", "", "", "")
weekxland_fe<-c(29, "weekxland_fe", "Yes", "", "", "")
adj_rsq<-c(30, "adj_rsq", round(summary(mo1_brid)$adj.r.squared, 3), "", "", "")
obs<-c(31, "obs", nobs(mo1_brid), "", "", "")

mode1_b<-rbind(mode1_b, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
mode1_b$model<-"Model1"
mode1_b$order<-as.numeric(mode1_b$order)
names(mode1_b)<-paste0(names(mode1_b), "_1")



mode2_b<-subset(model2, select = c("order", "term", "estimate", "std_er_star", "con_std_er_star"))
means<-c(23, "mean_mob", round(mean_sd_fun(mo8_brid)[1], 3), "", "", "")
sds<-c(24, "sd_mob", round(mean_sd_fun(mo8_brid)[2], 3), "", "", "")
no_cross_sections<-c(25, "cross_sections", round(week_obs_count(mo8_brid)[3], 3), "", "", "")
no_time_units<-c(26, "time_units", round(week_obs_count(mo8_brid)[1], 3), "", "", "")
county_fe<-c(27, "county_fe", "No", "", "", "")
week_fe<-c(28, "week_fe", "Yes", "", "", "")
weekxland_fe<-c(29, "weekxland_fe", "Yes", "", "", "")
adj_rsq<-c(30, "adj_rsq", round(summary(mo8_brid)$adj.r.squared, 3), "", "", "")
obs<-c(31, "obs", nobs(mo8_brid), "", "", "")

mode2_b<-rbind(mode2_b, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
mode2_b$model<-"Model2"
mode2_b$order<-as.numeric(mode2_b$order)
names(mode2_b)<-paste0(names(mode2_b), "_2")



mode3_b<-subset(model3, select = c("order", "term", "estimate", "std_er_star", "con_std_er_star"))
means<-c(23, "mean_mob", round(mean_sd_fun(mo1_brid_nokul)[1], 3), "", "", "")
sds<-c(24, "sd_mob", round(mean_sd_fun(mo1_brid_nokul)[2], 3), "", "", "")
no_cross_sections<-c(25, "cross_sections", round(week_obs_count(mo1_brid_nokul)[3], 3), "", "", "")
no_time_units<-c(26, "time_units", round(week_obs_count(mo1_brid_nokul)[1], 3), "", "", "")
county_fe<-c(27, "county_fe", "No", "", "", "")
week_fe<-c(28, "week_fe", "Yes", "", "", "")
weekxland_fe<-c(29, "weekxland_fe", "Yes", "", "", "")
adj_rsq<-c(30, "adj_rsq", round(summary(mo1_brid_nokul)$adj.r.squared, 3), "", "", "")
obs<-c(31, "obs", nobs(mo1_brid_nokul), "", "", "")

mode3_b<-rbind(mode3_b, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
mode3_b$model<-"Model3"
mode3_b$order<-as.numeric(mode3_b$order)
names(mode3_b)<-paste0(names(mode3_b), "_3")


mode4_b<-subset(model4, select = c("order", "term", "estimate", "std_er_star", "con_std_er_star"))
means<-c(23, "mean_mob", round(mean_sd_fun(mo8_brid_nokul)[1], 3), "", "", "")
sds<-c(24, "sd_mob", round(mean_sd_fun(mo8_brid_nokul)[2], 3), "", "", "")
no_cross_sections<-c(25, "cross_sections", round(week_obs_count(mo8_brid_nokul)[3], 3), "", "", "")
no_time_units<-c(26, "time_units", round(week_obs_count(mo8_brid_nokul)[1], 3), "", "", "")
county_fe<-c(27, "county_fe", "No", "", "", "")
week_fe<-c(28, "week_fe", "Yes", "", "", "")
weekxland_fe<-c(29, "weekxland_fe", "Yes", "", "", "")
adj_rsq<-c(30, "adj_rsq", round(summary(mo8_brid_nokul)$adj.r.squared, 3), "", "", "")
obs<-c(31, "obs", nobs(mo8_brid_nokul), "", "", "")

mode4_b<-rbind(mode4_b, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
mode4_b$model<-"Model4"
mode4_b$order<-as.numeric(mode4_b$order)
names(mode4_b)<-paste0(names(mode4_b), "_4")



mode5_b<-subset(model5, select = c("order", "term", "estimate", "std_er_star", "con_std_er_star"))
means<-c(23, "mean_mob", round(mean_sd_fun(mo1_bond)[1], 3), "", "", "")
sds<-c(24, "sd_mob", round(mean_sd_fun(mo1_bond)[2], 3), "", "", "")
no_cross_sections<-c(25, "cross_sections", round(week_obs_count(mo1_bond)[3], 3), "", "", "")
no_time_units<-c(26, "time_units", round(week_obs_count(mo1_bond)[1], 3), "", "", "")
county_fe<-c(27, "county_fe", "No", "", "", "")
week_fe<-c(28, "week_fe", "Yes", "", "", "")
weekxland_fe<-c(29, "weekxland_fe", "Yes", "", "", "")
adj_rsq<-c(30, "adj_rsq", round(summary(mo1_bond)$adj.r.squared, 3), "", "", "")
obs<-c(31, "obs", nobs(mo1_bond), "", "", "")

mode5_b<-rbind(mode5_b, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
mode5_b$model<-"Model5"
mode5_b$order<-as.numeric(mode5_b$order)
names(mode5_b)<-paste0(names(mode5_b), "_5")



mode6_b<-subset(model6, select = c("order", "term", "estimate", "std_er_star", "con_std_er_star"))
means<-c(23, "mean_mob", round(mean_sd_fun(mo8_bond)[1], 3), "", "", "")
sds<-c(24, "sd_mob", round(mean_sd_fun(mo8_bond)[2], 3), "", "", "")
no_cross_sections<-c(25, "cross_sections", round(week_obs_count(mo8_bond)[3], 3), "", "", "")
no_time_units<-c(26, "time_units", round(week_obs_count(mo8_bond)[1], 3), "", "", "")
county_fe<-c(27, "county_fe", "No", "", "", "")
week_fe<-c(28, "week_fe", "Yes", "", "", "")
weekxland_fe<-c(29, "weekxland_fe", "Yes", "", "", "")
adj_rsq<-c(30, "adj_rsq", round(summary(mo8_bond)$adj.r.squared, 3), "", "", "")
obs<-c(31, "obs", nobs(mo8_bond), "", "", "")

mode6_b<-rbind(mode6_b, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
mode6_b$model<-"Model6"
mode6_b$order<-as.numeric(mode6_b$order)
names(mode6_b)<-paste0(names(mode6_b), "_6")



mode7_b<-subset(model7, select = c("order", "term", "estimate", "std_er_star", "con_std_er_star"))
means<-c(23, "mean_mob", round(mean_sd_fun(mo1_bond_kul)[1], 3), "", "", "")
sds<-c(24, "sd_mob", round(mean_sd_fun(mo1_bond_kul)[2], 3), "", "", "")
no_cross_sections<-c(25, "cross_sections", round(week_obs_count(mo1_bond_kul)[3], 3), "", "", "")
no_time_units<-c(26, "time_units", round(week_obs_count(mo1_bond_kul)[1], 3), "", "", "")
county_fe<-c(27, "county_fe", "No", "", "", "")
week_fe<-c(28, "week_fe", "Yes", "", "", "")
weekxland_fe<-c(29, "weekxland_fe", "Yes", "", "", "")
adj_rsq<-c(30, "adj_rsq", round(summary(mo1_bond_kul)$adj.r.squared, 3), "", "", "")
obs<-c(31, "obs", nobs(mo1_bond_kul), "", "", "")

mode7_b<-rbind(mode7_b, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
mode7_b$model<-"Model7"
mode7_b$order<-as.numeric(mode7_b$order)
names(mode7_b)<-paste0(names(mode7_b), "_7")



mode8_b<-subset(model8, select = c("order", "term", "estimate", "std_er_star", "con_std_er_star"))
means<-c(23, "mean_mob", round(mean_sd_fun(mo8_bond_kul)[1], 3), "", "", "")
sds<-c(24, "sd_mob", round(mean_sd_fun(mo8_bond_kul)[2], 3), "", "", "")
no_cross_sections<-c(25, "cross_sections", round(week_obs_count(mo8_bond_kul)[3], 3), "", "", "")
no_time_units<-c(26, "time_units", round(week_obs_count(mo8_bond_kul)[1], 3), "", "", "")
county_fe<-c(27, "county_fe", "No", "", "", "")
week_fe<-c(28, "week_fe", "Yes", "", "", "")
weekxland_fe<-c(29, "weekxland_fe", "Yes", "", "", "")
adj_rsq<-c(30, "adj_rsq", round(summary(mo8_bond_kul)$adj.r.squared, 3), "", "", "")
obs<-c(31, "obs", nobs(mo8_bond_kul), "", "", "")

mode8_b<-rbind(mode8_b, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
mode8_b$model<-"Model8"
mode8_b$order<-as.numeric(mode8_b$order)
names(mode8_b)<-paste0(names(mode8_b), "_8")



#Step1: Collecting the models
list_order<-c(mode1_b$order_1, mode2_b$order_2, mode3_b$order_3, mode4_b$order_4,
              mode5_b$order_5, mode6_b$order_6, mode7_b$order_7, mode8_b$order_8)
#Step2: Ordering rows
list_order<-unique(sort(list_order))
list_order_df<-data.frame(list_order)
names(list_order_df)<-"order"
list_order_df$order<-as.numeric(list_order_df$order)


new_df2<-left_join(list_order_df, as.data.frame(mode1_b), by = c("order"="order_1"))

new_df3<-left_join(new_df2, as.data.frame(mode2_b), by = c("order"="order_2"))
new_df3$model<-new_df3$model_1
new_df3$model<-ifelse(is.na(new_df3$model), new_df3$model_2, new_df3$model)
new_df3$term<-new_df3$term_1
new_df3$term<-ifelse(is.na(new_df3$term_1), new_df3$term_2, new_df3$term)
new_df3<-subset(new_df3, select = -c(term_1, model_1, term_2, model_2))

new_df4<-left_join(new_df3, as.data.frame(mode3_b), by = c("order"="order_3"))
new_df4$model<-ifelse(is.na(new_df4$model), new_df4$model_3, new_df4$model)
new_df4$term<-ifelse(is.na(new_df4$term), new_df4$term_3, new_df4$term)
new_df4<-subset(new_df4, select = -c(term_3, model_3))

new_df5<-left_join(new_df4, as.data.frame(mode4_b), by = c("order"="order_4"))
new_df5$model<-ifelse(is.na(new_df5$model), new_df5$model_4, new_df5$model)
new_df5$term<-ifelse(is.na(new_df5$term), new_df5$term_4, new_df5$term)
new_df5<-subset(new_df5, select = -c(term_4, model_4))

new_df6<-left_join(new_df5, as.data.frame(mode5_b), by = c("order"="order_5"))
new_df6$model<-ifelse(is.na(new_df6$model), new_df6$model_5, new_df6$model)
new_df6$term<-ifelse(is.na(new_df6$term), new_df6$term_5, new_df6$term)
new_df6<-subset(new_df6, select = -c(term_5, model_5))

new_df7<-left_join(new_df6, as.data.frame(mode6_b), by = c("order"="order_6"))
new_df7$model<-ifelse(is.na(new_df7$model), new_df7$model_6, new_df7$model)
new_df7$term<-ifelse(is.na(new_df7$term), new_df7$term_6, new_df7$term)
new_df7<-subset(new_df7, select = -c(term_6, model_6))

new_df8<-left_join(new_df7, as.data.frame(mode7_b), by = c("order"="order_7"))
new_df8$model<-ifelse(is.na(new_df8$model), new_df8$model_7, new_df8$model)
new_df8$term<-ifelse(is.na(new_df8$term), new_df8$term_7, new_df8$term)
new_df8<-subset(new_df8, select = -c(term_7, model_7))

new_df9<-left_join(new_df8, as.data.frame(mode8_b), by = c("order"="order_8"))
new_df9$model<-ifelse(is.na(new_df9$model), new_df9$model_8, new_df9$model)
new_df9$term<-ifelse(is.na(new_df9$term), new_df9$term_8, new_df9$term)
new_df9<-subset(new_df9, select = -c(term_8, model_8))



new_df9$formal<-NA
new_df9$formal[new_df9$term=="bridging_tot_pop_dummy"]<-"Bridging"
new_df9$formal[new_df9$term=="bridging_nokul_tot_pop_dummy"]<-"Bridging (Excl. Culture)"
new_df9$formal[new_df9$term=="bonding_tot_pop_dummy"]<-"Bonding"
new_df9$formal[new_df9$term=="bonding_kul_tot_pop_dummy"]<-"Bonding (Incl. Culture)"
new_df9$formal[new_df9$term=="post_treat"]<-"Post Week 11 Dummy"
new_df9$formal[new_df9$term=="post_treat_bridging_tot_pop_dummy"]<-"Post Week 11 x Bridging"
new_df9$formal[new_df9$term=="post_treat_bridging_nokul_tot_pop_dummy"]<-"Post Week 11 x Bridging (Excl. Culture)"
new_df9$formal[new_df9$term=="post_treat_bonding_tot_pop_dummy"]<-"Post Week 11 x Bonding"
new_df9$formal[new_df9$term=="post_treat_bonding_kul_tot_pop_dummy"]<-"Post Week 11 x Bonding (Incl. Culture)"
new_df9$formal[new_df9$term=="lag_covid_pc"]<-"Lag Covid per Cap."
new_df9$formal[new_df9$term=="pct_turnout"]<-"Pct. Turnout"
new_df9$formal[new_df9$term=="log_pop_total"]<-"Log Pop. Total"
new_df9$formal[new_df9$term=="log_gdp_per_cap"]<-"Log GDP/capita"
new_df9$formal[new_df9$term=="pop_den"]<-"Pop. Density"
new_df9$formal[new_df9$term=="pct_college"]<-"Pct. College"
new_df9$formal[new_df9$term=="pct_pop_above_60"]<-"Pct. Pop. above 60"
new_df9$formal[new_df9$term=="pct_pop_under_35"]<-"Pct. Pop. below 35"
new_df9$formal[new_df9$term=="pct_service_sec_narrow"]<-"Pct. Empl. Hospitality and Transport"
new_df9$formal[new_df9$term=="pct_manufacturing"]<-"Pct. Empl. Manufacturing"
new_df9$formal[new_df9$term=="east_ger"]<-"East Germany"
new_df9$formal[new_df9$term=="gender_ratio"]<-"Gender Ratio"
new_df9$formal[new_df9$term=="pct_students"]<-"Pct. Students"
#new_df9$formal[new_df9$term=="uni_pop"]<-"University-Population Ratio"
#new_df9$formal[new_df9$term=="parks_pop"]<-"Parks-Population Ratio"
new_df9$formal[new_df9$term=="mean_mob"]<-"Mean Mobility"
new_df9$formal[new_df9$term=="sd_mob"]<-"SD Mobility"
new_df9$formal[new_df9$term=="cross_sections"]<-"No. Cross-Sections"
new_df9$formal[new_df9$term=="time_units"]<-"No. Time Units"
new_df9$formal[new_df9$term=="county_fe"]<-"County FE"
new_df9$formal[new_df9$term=="week_fe"]<-"Week FE"
new_df9$formal[new_df9$term=="weekxland_fe"]<-"Week X Land FE"
new_df9$formal[new_df9$term=="adj_rsq"]<-"Adj. R sq."
new_df9$formal[new_df9$term=="obs"]<-"Observations"



new_df10<-subset(new_df9, select = c(formal, estimate_1, std_er_star_1, con_std_er_star_1,
                                    estimate_2, std_er_star_2, con_std_er_star_2,
                                    estimate_3, std_er_star_3, con_std_er_star_3,
                                    estimate_4, std_er_star_4, con_std_er_star_4,
                                    estimate_5, std_er_star_5, con_std_er_star_5,
                                    estimate_6, std_er_star_6, con_std_er_star_6,
                                    estimate_7, std_er_star_7, con_std_er_star_7,
                                    estimate_8, std_er_star_8, con_std_er_star_8))


#Step3 Order columns
#final1<-new_df8[ , order(names(new_df8))]
final1<-new_df10
new_row<-c("Variable", rep(c("Estimate", "Robust SE", "Conley SE"), 8))
final2<-rbind(new_row, final1)


#Step7 centering coefficients
object_latex<-xtable(final2, type = "latex", 
                     caption = "Results",
                     #Note there are X columns to be aligned below.
                     #This is because of the index column.
                     align = paste(rep("l", ncol(final2)+1), collapse = ""))


comment <- list(pos = list(0), command = NULL)
comment$pos[[1]] <- c(nrow(object_latex))
comment$command <- c(paste("\\bottomrule\n",
                           "\\multicolumn{", ncol(final2)-1, "}{l} {\\parbox[t]{60cm}\n",
                           "{\\footnotesize \\textit{Note}: 
This table shows the effect of networks on mobility in the post-lockdown period. 
The dependent variable is mobility as measured by the government. 
The interactions between the period after week 11 and different types of networks (all networks, bridging, bonding)
are the variables of interest. All networks, bridging or bonding take the value of 1 for counties with 
density of these individual networks above the mean. Similarly, the post-lockdown period is marked with 1.
The standard errors are clustered at a county level in round brackets and 
adjusted for spatial (all networks - 455 km; bridging - 538 km; bonding - 455 km)
and serial (10 weeks) correlation in the square brackets. Significant at ∗10\\%, ∗∗5\\%, and ∗∗∗1\\%.}} \\\\ \n", sep = ""))


#Step8 Printing latex tex
tableLines<-print(object_latex,
                  #The following line controls where the lines are drawn
                  hline.after=c(-1, 1, 23),
                  floating = FALSE,
                  file = "./paper/tables/table_A7.tex", 
                  #caption.placement = "top",
                  add.to.row = comment,
                  booktabs = TRUE,
                  include.colnames=F,
                  include.rownames=FALSE,
                  sanitize.text.function=function(x){x})

multicolumns <- "& \\\\multicolumn{3}{c}{Model 1} & \\\\multicolumn{3}{c}{Model 2} & 
\\\\multicolumn{3}{c}{Model 3} & \\\\multicolumn{3}{c}{Model 4} & \\\\multicolumn{3}{c}{Model 5}
& \\\\multicolumn{3}{c}{Model 6} & \\\\multicolumn{3}{c}{Model 7} & \\\\multicolumn{3}{c}{Model 8} \\\\\\\\"
clines<-"\\\\cmidrule(lr){2-4} \\\\cmidrule(lr){5-7} \\\\cmidrule(lr){8-10} \\\\cmidrule(lr){11-13}
\\\\cmidrule(lr){14-16} \\\\cmidrule(lr){17-19} \\\\cmidrule(lr){20-22} \\\\cmidrule(lr){23-25}\\\\\\\\"
tableLines <- sub ("\\\\toprule\\n", paste0 ("\\\\toprule\n", multicolumns, clines, "\n"), tableLines) ## booktabs = TRUE
#tableLines <- sub ("\\\\hline\\n",   paste0 ("\\\\hline\n",   multicolumns, "\n"), tableLines) ## booktabs = FALSE
writeLines(tableLines, con = "./paper/tables/table_A7.tex", )


